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Abstract 

We investigate spatio-temporal structures in sheared polymer systems by solving a 
time-dependent Ginzburg-Landau model in two dimensions, (i) In polymer solutions 
above the coexistence curve, crossover from linear to nonlinear regimes occurs with 
increasing the shear rate. In the nonlinear regime the solution behaves chaotically 
with large-amplitude composition fluctuations. A characteristic heterogeneity length 
is calculated in the nonlinear regime, (ii) We also study dynamics of shear-band 
structures in wormlike micellar solutions under the condition of fixed stress. The 
average shear rate exhibits large temporal fluctuations with occurrence of large 
disturbances in the spatial structures. 

Key words: Two-fluid model, Shear-induced concentration fluctuations, 
Stress-diffusion coupling, Shear-banding, Rheological chaos. 



1 Introduction 

In complex fluids such as polymer solutions, liquid crystals, or surfactant sys- 
tems, the rheological properties can be very singular and complex. Application 
of flow field to complex fluids often gives rise to dramatic structural changes, 
sometimes coupled to phase separation [1,2]. Even chaotic spatio-temporal 
structures have been observed far away from equilibrium. 

In this paper we will present numerical results on two sheared polymer sys- 
tems: (i) In polymer solutions above the coexistence curve, application of flow 
has been observed to induce strong composition heterogeneities [3,4,5,6,7,8,9,10] 
The origin is now ascribed to a dynamical coupling between stress and dif- 
fusion in asymmetric viscoelastic mixtures [11,12,13,14]. This is in marked 
contrast to shear-induced homogenization or mixing observed in near-critical 
fluids [15,16] and in polymer mixtures with small entanglement effects [17,18]. 
Because some numerical results have already been presented [19,20,21,22], we 



Preprint submitted to Elsevier Science 



2 February 2008 



will here focus our attention on crossover from linear to nonlinear regimes with 
increasing the shear rate and chaotic behavior in the nonlinear regime. Fur- 
thermore, viscoelastic phase separation below the coexistence curve in shear 
flow has recently been studied, where polymer-rich domains are fragmented 
like broken gels [20] . (ii) Next we will briefly examine shear effects in wormlike 
micellar solutions where the stress-strain relation is nonmonotonic. Wormlike 
micells have significant chain flexibility and can break and combine [23]. In 
such systems shear bands composed of low-shear regions and high-shear re- 
gions parallel to flow have been observed [24,25,26,27] and Theological chaos, 
which occurs at low Reynolds numbers, has been discussed experimentally 
[28,29,30,31,32] and theoretically [33,34,35,36]. In the presence of shear bands 
macroscopic observables, such as the shear stress at fixed shear rate or the 
shear rate at fixed shear stress, exhibit irregular nonstationary chaotic behav- 
ior, which arises from large-scale disturbances of the shear bands. However, 
not enough understanding has yet been attained on such nonstationary behav- 
ior of the shear bands. In this paper we use the word chaotic loosely without 
rigorous characterization. 

Our numerical simulations will be performed in a Ginzburg-Landau scheme, 
where a conformation tensor representing the chain deformations is a dynamic 
variable. In previous papers numerical simulations (still in two dimensions) in 
this line have provided useful insights into the complex nonlinear shear effects 
[19,20,21,22]. The organization of this paper is as follows. We will explain our 
theoretical scheme in Section 2. We will then present numerical results for 
polymer solutions in Section 3 and for wormlike micellar solutions in Section 
4. 



2 Model equations in dimensionless forms 

Phase behavior of polymer solutions is usually described in terms of the poly- 
mer volume fraction only [37]. However, to describe the viscoelastic effects, 

it is convenient to introduce a symmetric tensor W = {W^} which represents 
chain conformations undergoing deformations. The free energy density / then 
depends on <p and W- For simplicity, we will measure the free energy density 
and the stress in units of a = k^T '/ 't> iV 3 / 2 where t> is the monomer volume 
and N is the monomer number on a chain. In terms of the scaled polymer 
density $ = iV 1//2 0, / reads 

/ = $ l n $ _ + 1^3 + || V$ |2 + 1&Q(W), (1) 

in the semidilute region <fi <C 1 close to the coexistence curve (in theta solvent). 
The space integral F = J drf gives the free energy. The first three terms 
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constitute the Flory-Huggins free energy and the fourth term represents the 
gradient free energy, where u is related to the interaction parameter \ by 
u = N l l 2 (2x — 1) [37]. We have u — 2 and $ = 1 at the critical point and u = 
$ + on the spinodal curve. We measure space in units of £ = bN 1 / 2 /2\/l8 
(~ gyration radius), where b = v^ d is the monomer length with d being the 
space dimensionality. The equilibrium average of Wij is the unit tensor 5ij and 
the deviation 5Wij = — 5ij gives rises to an increase in the free energy. We 

adopt the Gaussian form Q(W) = J2ij{^^ij) 2 ■ The coefficient g is of order 1 
in theta solvent and will be set equal to 1. 

Next we derive the dynamic equations using a two-fluid model. We assume that 
the polymer and the solvent have different velocities, v p and v s , respectively. 
The continuity equation is written as dcp/dt = — V • (</>v v ). In the semidilute 
region it is rewritten as 

4-$ + V- = -V- ($«;), (2) 
at 



in terms of the volume-averaged velocity v = 0v p + (1 — 0)v s and the relative 
velocity w = v p — v s . For <p <^ 1 we have v s = v and v p = v + w. The right 
hand side of Eq.(2) represents the diffusion. We assume that the solution is 

incompressible or V • v = 0. The motion of the conformation tensor W is 
driven by the polymer velocity v p as [1,22,38,39,40,41] 



(It + v p ■ - U(S ik 5W kj + 5W lk S kj ) + (sw ik n kj - n ik sw kj ) 
at k L 

= (Vj-upi + V iVpj ) - (1/r - VW^SWij, (3) 

where Vj = d/dxj, Sij = {ViV pj + VjV pi )/2, ft^ = {ViV pj - WjV pi )/2, and 
r is the very long rheological time. The left hand side is the frame-invariant 
time derivative of the tensor 5Wij and \a\ <_1 will be assumed in Section 4. 
The V is a kinetic coefficient originally introduced to explain shear banding in 
the Johnson-Segalman model [39,40], but is not indispensable in the presence 
of the gradient term in the free energy Eq.(l) for the composition. In our 
previous papers we have assumed a = 1 and V = for polymer solutions. 
However, we will assume a = 0.4 and V = 1 for wormlike micellar solutions. 
The network (viscoelastic) stress due to 5Wij consistent with Eq.(3) is given 
by [2] 



°pij = G E sw ik(hj + a5W kj ) + jGQSij, (4) 
k 4 
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where G = g$ 3 has the meaning of the shear modulus (in units of gq) for 
small rapid deformations. The linear viscosity is given by rj = i] + Gt in 
homogeneous states under weak steady shear, which is taken to be much larger 
than the solvent viscosity rj . Note that if G 3> 1, large viscoelastic effects 
are produced even for small SWij. The eigenvalues of SW^ remain small (at 
most of order 0.1) in our simulations [20], justifying the Gaussian form of Q. 
Physically, this means that the chains are weakly stretched in shear flow. In 
our case we simply have a pi j = GSWij from Eq.(4). 

Hereafter time will be measured in units of r = ?] /4o- (~ Qirr) £ 3 / k B T) . The 
velocities will be measured in units of £/tq = 4o"o^/^o- If we neglect the inertial 
effects (accelerations) in the two-fluid model, the velocities w and v read [20] 

w = -V^ + • a p - ^g$V$, v=(v)- ^[*w] ± , (5) 

in the semidilute region. The (v) is the mean shear flow, [• • -]± denotes taking 
the transverse part (whose Fourier component is perpendicular to the wave 
vector) and the inverse operation 1/V 2 may be expressed in terms of the 
Oseen tensor in the limit of large system size [2]. The stress force V • o" p in 
the diffusion flux &w gives rise to the so-called stress-diffusion coupling [2,14]. 
This means that diffusion is induced by stress imbalance V • o v ^ 0, which 
leads to various viscoelastic effects such as nonexponential decay in dynamic 
light scattering or shear-induced phase separation. 

We now have a closed set of dynamic equations for <3> and W ■ In particular, 
use of these equations yields the time derivative of the free energy F, 

If (v) = 0, we have dF/dt <0. The steady condition dF/dt = is realized in 
equilibrium where 5F/5& = const, and SW^ = 0. This is a self-consistent con- 
dition of dynamic equations which have stable equilibrium solutions [2] . Under 
this condition we may add the Langevin thermal noise terms in the dynamic 
equations Eqs.(2) and (3), although they have been omitted for simplicity. 

Assuming that the flow direction is along the x axis and the velocity-gradient 
direction is along the y axis, we give the expressions for the shear stress a xy 
and the normal stress difference Ni = a xx — a yy in units of ctq [20] , 

cr xy = <T pxy - 4$- 1 V :c <!>V ? ,<I> + A(V x v y + V y v x ), (7) 

iVl = iVpl - A<$>-\V X <&V X <& - Vy<S>Vy<S>) + 8(V X V X - VyVy), (8) 
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with v in units of £/tq. The first terms are the network stress components, 
the second ones arises from the gradient free energy (yielding the surface 
tension contributions in two-phase states), and the third ones are the viscous 
contributions. In homogeneous states with weak shear 7 = (V y v x ) we find 
& X y = (g& 3 T + 4)7 with r in units of r , so r]/r] = 1 + g$ 3 r/4 3> 1. For our 
inhomogeneous flows below, we also have (cr^) = (opr/) and (Ni) = (N pl ). 



3 Nonlinear structures in semidilute polymer solutions 

3. 1 Numerical method for sheared fluids 

For polymer solutions we set a — 1, T> — 0, and r = tq(0.2 + $ 4 )/(l + Q), 
where Q in the denominator gives rise to shear-thinning for 7T >1 even in 
homogeneous states [2,20]. For weak chain-stretching in the semidilute region, 
we have r/r = $ 4 ^> 1. Here we explain our numerical method. We integrate 
the dynamic equations in two dimensions on a 256 x 256 lattice. Our numerical 
scheme is to use the deformed coordinates, 

x' = x - nf(t)y, y' = y, (9) 

under the periodic boundary condition f(x', y' + L) = f(x' + L, y') = f(x', y') 
for any quantity f(x',y') in the deformed-space representation [42]. Then we 
can use the FFT method in the deformed space to calculate v from Eq.(5). 
Here we take 7(t) = it at constant 7. In Section 4, however, 7(2) will be 
determined such that the average shear stress (cr xy ) is held at a given constant. 
As the initial condition at t = we assign Gaussian random compositions at 
each lattice point, with mean value (<&) and variance 0.1. For t > we solve 
the dynamic equations without the Langevin noise terms. 



3.2 Mechanism of shear-induced fluctuation enhancement 

For very small shear rate 7 C r" 1 the linear stability analysis predicts the 
composition fluctuation enhancement above the coexistence curve owing to the 
stress-diffusion coupling [11,12,13,14]. At very long wavelengths the relaxation 
rate strongly depends on the direction of the wave vector k = k~ l k as [2] 



(10) 
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where D is the mutual diffusion constant, K os = o"o($— m$ 2 +$ 3 ) is the osmotic 
modulus, and ^ e (~ (v/VoY^b/^) is the viscoelastic length much longer than 
the correlation length (~ 6/0). The reduction of the relaxation rate (growth 
rate in spinodal decomposition) or the Onsager kinetic coefficient by the factor 
l/(l+/c 2 £ 2 e ) has been detected by scattering experiments without shear [43,44]. 
Note that T(k) < in the direction of k x = k y for K < 70(9r//0. In the 
semidilute region d In r]/d In ~ 6 experimentally (7 in our simulation since 
7] = Gt). It then follows the condition of strong fluctuation enhancement, 

7 > K os /6 V ~ O.lr- 1 , (11) 

in fair agreement with a scattering experiment [6]. Under Eq.(lO) the linear 
theory cannot be used and we encounter a nonlinear regime. The orign of 
the fluctuation enhancement is as follows: when the viscoelastic properties 
strongly depend on 0, appearance of inhomogeneities can lower the network 
free energy (oc g in Eq.(l)) and the stress (at fixed shear rate). By the same 
reason, when long and short polymers pass through a capillary tube, the longer 
ones migrate toward the tube center [14]. 

3. 3 Shear-induced fluctuations above the coexistence curve 



Figure 1 shows time evolution of the variance y((5$) 2 )(t), the shear stress 
(o~ xy }(t), and the normal stress difference (Ni)(t) after application of shear 
at t = for 7 = 0.05 = 0.8/r. The averages are taken over all the lat- 
tice points. Here the system is above the coexistence curve with u — 2 and 
($) = 2. Though there is no thermal noise in the dynamic equations, the 
variance saturates at about 0.6 in the late stage. Figure 2 suggests that the 
system is in a state of deterministic chaos, undergoing partial or incomplete 
phase separation. Both (cr xy )(t) and (Ni)(t) exhibit a pronounced peak at 
the onset of fluctuation enhancement at t ~ 25 and noisy behavior follows in 
the subsequent dynamical steady state. The noise of these quantities arises 
from breakup of the polymer-rich domains supporting the stress. The effective 
viscosity r] e s = {cxy)/i m the steady state is about 8 times larger than the 
solvent viscosity i]q. It would be 28?7o for the homogeneous case, so marked 
shear thinning is realized. For comparison the shear stress (upper broken line) 
and the normal stress difference (lower broken line) are calculated for the 
homogeneous case. 

Figure 2 displays the normalized composition <&(r,t) at t = 45 (left) and 
t = 300 (right) in the run which produced Fig.l. In the top figures we cannot 
see well-defined interfaces, but the composition changes over a mesoscopic, 
characteristic length in the steady state. We numerically find that the 
average gradient term (2|V$| 2 /$) and the average elastic energy (g<fr 3 Q/4) in 



6 



«t»=2 u=2 g=1 Y=0.05 



Homogeneous 




300 



Fig. 1. Time evolution of yj (5& 2 ) (dotted line), {cr xy ), and {N±) (solid lines) af- 
ter application of shear. At t = 0, the system is above the coexistence curve. The 
shear-induced fluctuations (shown in Fig.2) give rise to the pronounced stress over- 
shoots. The upper curves are those for the homogeneous case. 

the dimensionless free energy density Eq.(l) are of the same order for various 
7 and ($). We estimate the former as ( and the latter as (a xy ) 2 /($) 3 
to obtain 



R h ~($) 2 /(a xy ), (12) 

or i?h ~ bG/(f>(a xy ) in the original units. The behavior of (<J xy ) as a function 
of 7 can be known from Fig. 5. On the other hand, to the stress, the nework 
contributions are mostly much larger than the gradent contributions in Eqs.(7) 




x t=45 t=300 



Fig. 2. <&(x,y,t) at t = 45 (left) and t = 300 (right) above the coexistence curve. 
The top figures are the profiles along the x axis at y = 128. The parameters are the 
same as those in Fig.l. For t > 100 a dynamical steady state is realized. 
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and (8). 



The length should be measurable by scattering. In Fig. 3 the time- averaged 
structure factor S(k x ,k y ) = (|$k| 2 ) in the dynamical steady state is shown. 
It has sharp peaks in the k x -k y plane close to the k x axis, presumably at 
k x ~ ±27r/i?h- This is in agreement with the experimental structure factor for 
7r ~ 1 in the k x -k y plane [3]. It is worth noting that the light scattering in the 
k x -k z plane was also intensified along the /c x -axis for jr ~ 1 [6]. Microscope 
pictures in the x-z plane of the shear-induced composition fluctuations closely 
resemble the pattern at t — 300 in Fig. 2 [8]. These suggest stronger correlations 
perpendicular to the flow direction than along the flow for 7T ~ 1. Below the 
coexistence curve also, the structure factor exhibits similar behavior [20]. 



Fig. 3. Bird-view (left) and contour plot (right) of the time-averaged structure 
factor S(k x ,k y ) = (|$k| 2 ) in the dynamical steady state in Fig.2. 



3.4 Chaotic behavior above the coexistence curve 

As shown in Fig. 4, there emerge two dynamic regimes at longer times. When 
7 is smaller than a characteristic shear rate 7 n i(— 0.01 for (<3>) = 2), the fluc- 
tuations enhanced in the early stage slowly relax to zero due to the convective 
deformation by shear in the absence of the thermal noise. In this linear or ther- 
mal regime, the steady-state fluctuations can be attained only in the presence 
of the thermal noise. It should be noted that in the previous theories the fluc- 
tuation intensity was calculated in the linear approximation in the presence 
of the thermal noise. [11,12,13,14]. In the nonlinear or chaotic regime 7 > 7^, 
the large-amplitude fluctuations are continuously produced by the nonlinear 
interactions among them even without the thermal noise. In Fig.5 we show a 
discontinuous change in the steady-state variance with varying 7 for three (<&) 
(left), which occurs roughly at 7 ~ 0.2/r with r = r (0.2 + $ 4 ). The variance 
goes to zero in the linear regime in the absence of the thermal noise. Of course, 
the changeover becomes continuous in the presence of the thermal noise. In 
Fig.5 we also show the normalized effective viscosity Ves/vo — foxy)/ Vol versus 
7 (right), which exhibits Newtonian behavior in the linear regime and strong 
shear-thinning behavior with i] c g ~ I/7 in the nonlinear regime. 
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Fig. 4. Time evolution of ((5<&) 2 )(t) for four values of 7 without the Langevin 
noise. A chaotic steady state is realized for the upper three curves. For the smallest 
shear 7 = 0.01 the fluctuations enhanced in the early stage slowly relax to zero. 
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Fig. 5. Steady-state variance \J ((5<&) 2 ) (left) and normalized effective viscosity 
ijes/vo (right) versus 7 for three values of ($). Here there is no Langevin noise 
terms and a sharp discontinuous transition occurs between the linear and nonlinear 
regimes. Shear-thinning behavior is marked in the nonlinear regime. 

We give some remarks, (i) In one-phase states of near-critical fluids without 
viscoelasticity, the critical fluctuations are always suppressed by shear flow 
below the equilibrium level and the thermal noise is indispensable to produce 
the steady-state structure factor [2,15]. (ii) Below the coexistence curve, on 
the contrary, the random noise terms are generally irrelevant in producing 
large-scale phase-separation patterns and are usually omitted in simulations 
without and with shear flow [2] . However small the shear rate is, the domains 
are eventually deformed nonlinearly by applied shear flow as the domains 
grow. Therefore, the linear shear regime is nonexistent in late-stage phase 
separation. For Newtonian fluid mixtures with small surface tension, string- 
like domain structures can be formed in steady states [18]. (iii)The system 
considered here has very low Reynolds numbers, so that the fluid inertia plays 
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no role in the irregular patterns. 



4 Nonlinear structures in wormlike micellar solutions 



4-1 Nonmonotonic stress-strain relation 




Fig. 6. Stress-strain relation Eq.(12) for wormlike micellar solutions in homogeneous 
states. The region d/Jxy/dj < is unstable and shear bands appear. 

We present preliminary results on sheared wormlike micellar solutions where 
a = 0.4, V — 1, and r = r (0.2 + $ 4 ) in Eq.(3). Several authors have used 
the framework of the two- fluid dynamics including the composition [22,41]. For 
\a\ < 1 it is important that the stress-strain curve is nonmonotonic. That is, for 
homogeneous shear flow, Eq.(3) is solved to give 5W xy = 7t/[1 + ( A fr) 2 (l — a 2 )} 
and 5W xx /(a + 1) = SW yy /(a - 1) = irb~W xy . Then Eq.(4) gives 



. ^. [l + ( 7 r) 2 (l + a 2 )] 

o xy = Vol + G 7 r V ( ' i^r (13) 

[l + (7r) 2 (l-a 2 )] 

For rjo/Gr ^ 1 the curve has a maximum at 7 ~ 1/r and a minimum at 
7 ~ (G/rjor) 1 / 2 for a not very close to 1. Fig. 6 displays the curve for Gr/r]o = 
32.4. It is known that homogeneous states are unstable for da xy jd^i < 0. Shear 
banding instability occurs in the region of a xy where 7 is multi-valued. 



4-2 Numerical method at fixed stress 



For this system there have been experiments both at constant shear rate 
[28,29,30,31] and at constant shear stress [32]. Since simulations of a model 
similar to ours have already been performed at constant shear rate for a = 
[22], we here show results at constant shear stress. We solve the dynamic 
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equations in Section 2 on a 128 x 128 lattice using the numerical scheme in 
Section 3 [42]. We impose the periodic boundary condition in the deformed 
space, so the boundary effects do not appear in our scheme. We use the Stokes 
approximation V • a — to find / L dxV y a xy (r, t) = 0, where a is the (total) 
stress tensor and L is the system length in the x axis. Then the lateral average 
defined by 

l 

E = y dx-a xy (r,t) (14) 
o 



becomes independent of y. If E is a constant independent of time, Eq.(7) yields 
the average shear rate dependent on time in the dimensionless form, 



4<7>(i) = E - (a pxy - 4f- 1 V I $V J $>, (15) 

where ( ) denotes taking the space average. At each time step we calculate 
the average strain fo dt' (i) (t') as the shear strain ^(t) in Eq.(9). 



4-3 Chaotic behavior with concentration inhomogeneities at fixed stress 




Fig. 7. Time evolution of the average shear rate (7)(i) for (a xy ) = 5.4 and ($) = 2.0. 

Because we set u = —4, the system is in a one-phase state before application of 
shear. Figure 7 presents the time evolution of the average shear rate (7) (t) = 
(dv x /dy) in the initial stage of band formation after application of shear flow 
at t — 0. In Fig. 8 we show snapshots of <&(r,£) and 7(7*, t) = V y v x at various 
times in this stage. We can see that the homogeneous state becomes unstable at 
t ~ 400 and the shear bands subsequently begin to grow in the time region t = 
400 ~ 500. The pseudo one-dimensional band structures are soon deformed in 
the flow direction, resulting in the irregular oscillatory behavior of {j)(t). 
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t=400 450 500 600 




Fig. 8. Snapshots of $(r,t) and i(r,t) = V y v x at t = 400, 450, 500, 600. The dark 
regions take larger values than the bright regions. 

We next explain the long-time behavior. Figure 9 shows the chaotic temporal 
fluctuations of the space average (7)(t) = (dv x /dy) for (a xy ) = 5.4 and (<&) = 
2.0. We notice that frequently exhibits large amplitude fluctuations, 

which arises from large-scale disturbances of the interfaces between different 
shear bands. Figure 10 illustrates $(r,t) and 7(r*,i) = V y v x at t = 9000 
(A), t — 14000 (B), and t = 15250 (C) in Fig.9. The bird views of (A) and 
(C) are given also. We find that the system jumps between the two states, 
(A) and (B), alternatively. In these states the thickness of the low-shear band 
(white) is considerably different. Furthermore, at t ~ 15000 in (C) the interface 
between two bands becomes globally unstable, resulting in abrupt expansion 
of the low-shear band. We can see that the velocity and the composition are 
strongly coupled. Thus these chaotic spatio-temporal behavior arise from the 
stress-diffusion coupling in our model. 



<$>=2 u=-4 g=1 <a xy >=5.4 




Fig. 9. Time evolution of the average shear rate (7)(t) for (a xy ) = 5.4 and ($) = 2.0. 
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(A) (B) (C) 




Fig. 10. Snapshots and bird views of $(r,t) and j(r,i) = V y v x at t = 9000 (A), 
t = 14000 (B), and t = 15250 (C) in Fig.9. The dark regions take larger values than 
the bright regions in the upper snapshots. 

5 Summary 

We have investigated complex spatio-temporal structures in sheared polymer 
systems by solving a time-dependent Ginzburg-Landau model. Though per- 
formed in two dimensions, we have revealed a number of unique features of 
sheared polymer systems. We summarize our main results, 
(i) In polymer solutions above the coexistence curve, Fig. 4 shows two dynamic 
regimes. In the linear regime 7 < 7 n i the shear-induced fluctuations can be 
maintained only in the presence of the thermal noise, while in the nonlin- 
ear regime 7 > 7 n i chaotic dynamical steady states are realized due to the 
nonlinear interactions among the fluctuations. The previous linear calcula- 
tions [11,12,13,14] are meaningful only in the linear regime, and our results in 
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Figs. 1-5 are those in the nonlinear regime. The stress overshoot in Fig.l occurs 
due to appearance of the shear-induced composition fluctuations in Fig.2. In 
the dynamical steady states of the nonlinear regime the heterogeneous devi- 
ation 5(j) is even on the order of the average ((f)) varying on the scales of the 
characteristic size ~ Rh as in the top figures of Fig.2. We have related Rh to 
the average shear stress as in Eq.(12) in the dynamical steady states, 
(ii) We have also presented our preliminary results of shear-band structures 
in wormlike micellar solutions under the condition of fixed stress. As shown 
in Figs. 9 and 10, the average shear rate exhibits large temporal fluctuations. 
However, we have not yet understood the mechanism why the shear bands 
exhibit irregular behavior in shear flow. 

We will shortly report simulations of wormlike micellar solutions both at fixed 
shear rate and at fixed shear stress. We are not aware of any previous simu- 
lations at fixed shear stress. As other problems not treated here, we mention 
rheological turbulence at low Reynolds numbers in dilute polymer solutions 
[45] and defect turbulence in nematic liquid crystalline polymers in shear flow 
[2,46,47]. Both shear banding and tumbling have been observed in suspensions 
of rods [48]. Liquid crystalline molecules can tumble in shear flow, constitut- 
ing a coupled oscillator system [49], and loss of their phase coherency leads 
to proliferation of disclinations. Thus a number of nonequilibrium problems 
remain largely unexplored in complex fluids, which could also be interesting 
subjects in nonlinear science and in fluid mechanics. 
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